Experimental and numerical study on gas production decline trend under ultralong-production-cycle from shale gas wells

It is of engineering interest to explore recovered shale gas composition and its effects on total gas production trend over a long-term extraction period. However, there are previous experimental studies mostly focused on short term development for small scaled cores, which is less convincing to mimic reservoir-scaled shale production process. In addition, the previous production models mostly failed to account for comprehensive gas nonlinear effects. As a result, in this paper, to illustrate the full-life-cycle production decline phenomenon for shale gas reservoir, dynamic physical simulation was performed for more than 3433 days to simulate shale gas transport out of the formations over a relatively long production period. Moreover, a five-region seepage mathematical model was then developed and was subsequently validated by the experimental results and shale well production data. Our findings show that for physical simulation, both the pressure and production declined steadily at an annual rate of less than 5%, and 67% of the total gas in the core was recovered. These test data supported earlier finding that shale gas is of low flow ability and slow pressure decline in the shale matrices. The production model indicated that free gas accounts for the majority of recovered shale gas at the initial stage. Based on a shale gas well example, free gas extraction makes up 90% of produced total gas. The adsorbed gas constitutes a primary gas source during the later stage. Adsorbed gas contributes more than 50% of the gas produced in the seventh year. The 20-year-cumulative adsorbed gas makes up 21% of the EUR for a single shale gas well. The results of this study can provide a reference for optimizing production systems and adjusting development techniques for shale gas wells throughout the combinations of mathematical modeling and experimental approaches.

China is renowned for its abundant shale gas resources, with total recoverable reserves estimated at 21.84 trillion m 3 . In the southern region of China, recoverable marine facies shale gas resources amount to 14.1 trillion m 3 . Nearly a decade of exploration and development efforts has enabled China to commercialize shale gas development. In 2021, China's total shale gas production reached 22.8 billion m 31 . Integrated geological and engineering technologies for reservoirs shallower than 3500 m have gradually taken shape, while the primary technologies used to explore and develop shale gas (e.g., horizontal well drilling, staged volume fracturing, and factory-like operation technologies) have basically been established. Nevertheless, there are still some urgent problems to be solved on development technologies and seepage theories for shale gas systems. For example, shale gas production varies appreciably between wells and declines in a yet-to-be-established pattern, while the available shale gas development systems require optimization. These problems, to some extent, hinder the efficient, sustained development of shale gas 2 .
Shale gas reservoirs differ from conventional gas reservoirs in that shale function as both the source rocks that generate natural gas and the reservoir and cap rocks that host and preserve natural gas. The Sichuan Basin and its peripheral areas-China's current main shale gas-producing region-provide a marine facies deep-water continental-shelf sedimentary environment. Here, high-quality shale reservoirs distributed continuously over a large area have developed in the Lower Silurian Longmaxi Formation (LMX Fm). Characteristically, the shale www.nature.com/scientificreports/ proportions during production period. The results of this study can provide important theoretical and practical guidance for efficiently developing shale gas wells.

Physical simulation experiment method
A full-diameter core sample of the Longyi 1 sublayer of the LMX Formation retrieved from well Y10 in Zhaotong in southern Sichuan, with the collection depth of 2450 m, was used in this test. The length and diameter measurements of the core were 15 cm and 10 cm, respectively. To maximally preserve its original characteristics, the core was quickly cleaned and then placed in a pressure-retaining coring device immediately after retrieval from the borehole while on site, followed by saturation with natural gas to 10 MPa. The valves at both ends of the core were then closed before it was transported to the laboratory for analysis. A parallel sample design was used in the test. Samples retrieved from the areas adjacent to the area where the full-diameter core was collected were tested to determine the physical properties of the shale and were also used for comparison (see Table 1 for the results). For example, some samples were ground to powder that could pass through sieves with mesh sizes of 60-100. The powder was subsequently subjected to an isothermal adsorption test to determine the adsorbed gas content. Because the shale in the study area contains developed beds, core plugs were prepared using the wirecutting technique to inhibit microfracture initiation during the conventional string drilling process. Methane (purity: 99.99%) was used in the test to simulate the shale gas flow in the reservoir. A GAI-100 high-pressure gas adsorption isotherm system (Core Laboratories, the United States) was used to measure the adsorbed gas content. The helium expansion method was adopted to measure the porosity on a PoroPerm-200 pore permeability measurement system. The classical steady-state permeability measurement method was employed in tandem with Klinkenberg's theory (included to account for the slip effect) to calculate the corresponding Klinkenberg permeability. Table 1 summarizes the basic parameters of the core used in the test. Evidently, both the porosity and permeability of the shale in the reservoir are very low, averaging approximately 2% and 0.00005 mD, respectively. With a flow conductivity on the order of 100 mD m (considerably higher than the permeability of the matrix), artificial fractures in a shale gas reservoir can be regarded as infinite conduction fractures. Therefore, the gas provided by the matrix can be considered to flow in a one-dimensional (1D) manner from the matrix to the nearby fractures (Fig. 2). A method was established to physically simulate 1D gas flow based on the core collected in this study from the reservoir. This method was subsequently used to examine the capacity of the matrix to provide gas to the fractures at the formation pressure on a test setup developed in-house to physically simulate the coupled desorption, diffusion, and seepage of shale gas (Fig. 3). The main components of the setup included a Teledyne ISCO pump, a high-pressure core holder, intermediate containers, a confining pressure pump, a gas  www.nature.com/scientificreports/ mass flowmeter, and high-precision pressure sensors. Figure 3 shows the test procedure described in detail as follows. (1) The core holder containing the core retrieved from the field was connected to the simulation test setup. The high-pressure pump was used to increase the natural gas pressure in the high-pressure intermediate containers to 28 MPa and maintain it at this level. The core was saturated again and pressurized to the formation pressure level (28 MPa) through injection of methane at both ends, while the confining pressure was simultaneously increased to 50 MPa. As a corollary of the extremely high density and adsorption behavior of the shale, fully restoring the core to its in situ state required a very long time. The pressurization and saturation process continued for 200 days and was terminated when the pressures in the intermediate containers at the inlet and outlet of the core holder no longer decreased. Under this condition, the shale gas was almost completely restored to its original state in the host formation. Subsequently, the inlet and outlet control valves (valves 1 and 2, respectively) were closed, followed by the removal of the high-pressure gas source. Thus, the setup was ready for a full-life-cycle production simulation test. (2) Valve 2 was opened to commence the test. To prevent an excess gas flow rate at the initial stage and maintain a reasonable back pressure at the outlet, the back pressure was adjusted to the atmospheric pressure when the flow rate fell within a suitable range. Under this condition, full-life-cycle long-term gas production from the shale gas well was dynamically simulated. (3) Throughout the simulation test, the inlet and outlet pressures and the outlet flow rate were continuously measured in real time using the inlet and outlet pressure sensors (sensors 1 and 2, respectively) and the outlet gas mass flowmeter, respectively, and the data were collected using a computer data acquisition system. To improve the test accuracy, the gas flow was measured using the inlet gas mass flowmeter during the initial stage when the flow rate was relatively high, combined with use of the bubble or micropipe drainage method during the low-production stage.

Physical simulation test results and discussion
Isothermal adsorption results and pore structure characterization. Figure 4 shows the isothermal adsorption curve experimentally obtained by testing dry samples using the adsorption isotherm system. As is presented in the Fig. 3, the excess adsorption curve is obtained based on the high pressure absorption characterization model 18 , written in the Eq. (1): where V a is excess gas adsorption volume m 3 /kg; V L is Langmuir volume, m 3 /kg; P L is Langmuir pressure, Pa; β is modification factor caused by high pressure and high temperature adsorption, m 3 /kg; P is gas pressure in seepage field, Pa; ρ g is gas density, kg/m 3 . Observation of the measured isothermal adsorption curve in Fig. 4 shows the following. The adsorption (i.e., excess adsorption) increased sharply as the pressure increased before reaching 10 MPa. As the pressure increased from 10 to 15 MPa, the excess adsorption increased slowly (a decrease in the excess adsorption was observed in some samples). Further increasing the pressure beyond 15 MPa reduced the excess adsorption, as observed in all the samples. This finding is consistent with that obtained from a high-pressure supercritical isothermal adsorption test. The typical pattern of the high-pressure supercritical isothermal adsorption curve of the shale differs from the monotonic increase pattern of its low-pressure isothermal adsorption curve, suggesting that the excess adsorption was measured, instead of the absolute adsorption (i.e., the actual adsorption). An absolute adsorption correction model was used to correct the excess adsorption curve. An absolute adsorption curve (the red curve in Fig. 4) was thus obtained. The high-pressure isothermal adsorption pattern of the shale, characterized by an increase and a subsequent decrease in the excess adsorption with an increase in the pressure, is difficult to describe with the Langmuir model (used to describe absolute adsorption) and other subcritical models 29 . Therefore, a new model must be established. Observation of Eqs. (2) reveals the following relationship between excess and absolute adsorption: where G ex is the excess absorbed gas accumulated in the shale rock, m 3 ; G abs is amount of adsorbed gas accumulated in the shale formation, m 3 ; ρ g is free gas density, kg/m 3 ; ρ a is absorbed gas density, kg/m 3 .
An assumed adsorbed-phase density is required to use the equation presented above. Researchers have fitted excess adsorption curves with the liquid-phase density (423 kg/m 3 ), van der Waals density (373 kg/m 3 ), and critical density as the adsorbed-phase density 18,30 . In this study, a revised Langmuir model was used to calculate the adsorbed gas content.
Pore structure characteristics. A nitrogen adsorption analyzer was used to generate a nitrogen adsorption curve for the samples, as shown in Fig. 5. Evidently, the adsorption-desorption curve is not prominent. Observation of the hysteresis loop reveals that the adsorption increased rapidly as the relative pressure approached 1. The hysteresis loop also appears relatively small. These findings mainly reflect the presence of parallel plate pores with four open sides. Moreover, the shape of the hysteresis loop conforms to that of a Type H4 hysteresis loop, indicating the presence of slit pores. Relatively well-developed pores in various diameter ranges (from micropores to large pores) were all observed, and they were relatively well interconnected. Analysis of the test results for the samples reveals a specific surface area of 12.81 m 2 /g, a Barrett-Joyner-Halenda (BJH) total pore volume of 0.0276 mL/g, and an average pore diameter of 9.38 nm. www.nature.com/scientificreports/ Focused ion beam (FIB) is a technique used to achieve microbeam fabrication with an ion beam focused to a sub micrometer-or even nanometer-sized spot through a deflection system. Scanning electron microscopy (SEM) is employed to produce images of a sample. In this technique, various types of signals (e.g., secondary and backscattered electron signals) are first generated from the sample through the interactions between high-energy electrons and matter and are subsequently converted sequentially and proportionally through a detector into video signals that are then amplified and adjusted for point brilliance to form an image. In this study, a Helios NanoLab 650 dual beam SEM system (ion beam voltage: 500 V-30 kV; electron beam voltage: 350 V-30 kV) was used to produce two-dimensional microimages of the core samples, with the goal of revealing the pore structure (laminae and microfractures), nanoscale particles and pores, and mineral composition on the sample surfaces.
Inspection of the SEM images shows an extensive presence of pores in the OM. As seen in Fig. 6a, a large number of circular or oval-shaped bubble pores of a relatively uniform size were present in the kerogen. Most of these pores exhibited a spongy or honeycomb-like appearance with a diameter less than 100 nm. Figure 6b shows organic pores in a mixture of OM and clay. Changes in geological conditions during accumulation and thermal evolution lead to the development of multitudinous micropores or microcracks in OM. Distributed primarily along microbedding planes or sedimentary discontinuities, OM tends to form interconnected pore networks and thus has a relatively high permeability. Banded OM, pyrite, and a large amount of quartz and calcite are visible in Fig. 6c. The pores in this sample developed predominantly within the OM and were relatively poorly connected, conforming to the low permeability of the matrix.
Production pressure and production decline curves. Figure 7 shows the inlet pressure and gas production rate curves obtained from the test after 3356 days (September 15, 2012, through November 30, 2021).
Analysis of the production curve in Fig. 7 and Fig. 8 reveals the following. The tested core produced gas at a rate that was high but decreased rapidly during the initial stage and subsequently quickly stabilized then remained stable for a protracted period of time. The production decline curve closely resembles a single-well production decline curve. Gas was produced at a high rate during the first 40 days. Throughout this stage, the pressure decreased from 28 to 20 MPa, while the gas production rate decreased rapidly from 140 to less than 10 mL/d. Subsequently, there was a decrease in both the rates at which the daily gas production decreased and  www.nature.com/scientificreports/ the cumulative production increased. The pressure and daily gas production decreased to 18 MPa and 2 mL/d, respectively, after 200 days of production, followed by a long-lasting stage of low but stable production, during which the gas production rate remained stable and decreased extremely slowly. By the end of the production period examined in this study (i.e., the 3356th day), the inlet pressure and daily gas production had stabilized at 6.8 MPa and approximately 0.2 mL, respectively, while the cumulative gas production had reached 4837 mL. According to the calculation of the gas content, the core contained 6710 mL of gas in total. Taken together, these results indicate that 72.1% of the gas had been recovered. The extremely low permeability and adsorption and desorption behaviors were the primary causes of the slow pressure decline and a long period of low but stable gas production observed in the shale. Figure 9 shows a plot of the apparent pressure versus the cumulative gas production. Observation of Fig. 9 shows a relatively notable linear relationship between the cumulative gas production and apparent pressure with a decrease the apparent pressure (as indicated by the red dots in Fig. 9). The apparent formation pressure decreased further below 12 MPa at a decreasing rate, as indicated by the deviation from the blue straight line in Fig. 9, the extent of which increased as the apparent formation pressure decreased. This finding suggests a critical desorption pressure of 12 MPa for the adsorbed gas in the shale core at room temperature, which is consistent with that obtained from a high-pressure adsorption test on shale samples collected from the study area 16 . After the apparent formation pressure decreased below the critical desorption pressure, the adsorbed gas began to desorb extensively and contribute to the gas production, while the cumulative gas production surpassed the value calculated using the material balance equation. The difference in the cumulative gas production increased as the formation pressure decreased. In other words, as the formation pressure decreased, the desorbed gas content of the gas produced from the core increased. The cumulative gas production was also basically linearly related with the apparent formation pressure during the low and stable production stage.
According to the material balance equation, for a closed gas reservoir, if the adsorbed gas is not considered, then the apparent formation pressure decreases linearly as the free gas production increases, as shown in the equation below 35 :  www.nature.com/scientificreports/ where G 0 is total free gas reserves in the shale formations, m 3 ; G free is free gas amount at the reservoir pressure of P , m 3 ; Z i is gas compression factor at the original reservoir pressure of P i , dimensionless; Z is gas compression factor at the reservoir pressure of P , dimensionless; P i is initial pressure in the gas reservoir, Pa; P i is reservoir pressure at time t, Pa. Observation of the equation presented above shows that G 0 Zi/Pi and G 0 (i.e., total free gas content) are the slope and intercept of the cumulative shale gas production-apparent pressure curve for the early stage, respectively. Therefore, G 0 can be determined from the slope of the curve. The dotted red line in Fig. 9 signifies the calculated material balance curve for the free gas in the closed reservoir, which reflects the dynamic free gas production from the shale core and can be used as reference data to determine the free and adsorbed gas content of the gas produced at each formation pressure. The difference (approximately 1416.6 mL) between the free gas content calculated using Eq. (2) (4956 mL) and the gas content (i.e., the total free and adsorbed gas content) corresponding to the measuring point (6372 mL) is the adsorbed gas content of the gas produced from the core.
The free and adsorbed gas production patterns of the core can be further derived from the formation pressure and the patterns of the variation in the cumulative and daily gas production with the decrease in the formation pressure in conjunction with the basic theories of gas reservoir engineering 36,37 . The free gas content of the shale sample is obtained based on the pore volume and volume coefficient, while its adsorbed gas content is determined using a revised Langmuir adsorption Eq. 16 Let Pi and P denote the original formation pressure and the formation pressure at time t, respectively. Then, the amounts of free and adsorbed gas recovered from the core are calculated below: Figure 8. Measured cumulative production and gas production decline rate curves of the full-diameter core.

Figure 9.
Measured relationship between the cumulative gas production and P/z of the full-diameter core. www.nature.com/scientificreports/ where G free is free gas production at time t, m 3 ; G ab is adsorbed gas production at time t, m 3 ; V is shale reservoir volume, m 3 ; ∅ M is matrix porosity, dimensionless; Tsc is standard temperature, K; Psc is standard atmospheric pressure (0.101 MPa); T is reservoir temperature, K; ρ gi is gas density under initial reservoir conditions, kg/m 3 . Therefore, the cumulative amount of gas extracted from the core is Use of the model yields a total gas content of 6703 mL for the core, of which adsorbed and free gas account for 1644 (24.5%) and 5059 mL, respectively. A comparison with the gas content results derived from the fitted curve in Fig. 9 shows the following. Of the three values, the free gas content yielded by the model is the closest to the result derived from the curve, with a difference of only 2%. The total gas content yielded by the model differs relatively markedly (by 331 mL) from the result derived from the curve, primarily due to the difference in the adsorbed gas content (228 mL; approximately 13%). This disparity can be mainly ascribed to the error in the total gas content derived from the linearly fitted curve (black portion) in Fig. 8, which occurs because as the pressure decreases, adsorbed gas desorbs at an increasing rate, resulting in a further deviation of the black line from the line denoting the linear relationship and thus a total gas content greater than 6372 mL. Therefore, the amounts of adsorbed and free gas calculated based on the porosity and isothermal adsorption curve are reliable.
Inspection of Fig. 10 shows that the cumulative gas production curve basically overlaps the free gas production curve and is highly linearly related with the apparent mean pressure curve during the initial high-pressure stage. This phenomenon can be attributed to the nonsignificant impact of desorbed gas due to its relatively small quantity during this stage. As the pressure decreased, the adsorbed gas production increased and contributed more to the cumulative gas production, resulting in a deviation of the cumulative gas production curve from the straight line. In contrast, a high level of linearity between the free gas production and the apparent pressure persisted throughout the observation time. By the end of the production period examined in this study, the cumulative gas production had reached 4837 mL, of which free and adsorbed gas accounted for 3363 mL (71% of the total amount of free gas) and 1204 mL (73% of the total amount of adsorbed gas), respectively. The extraction of free gas was responsible for the high production during the initial stage, which was followed by the extraction of adsorbed gas. The physical simulation satisfactorily reflects the characteristics of shale gas production-high initial production, followed by a rapid decrease and stable production during the later stage. The contribution of adsorbed gas became pronounced only during the later stage, which led to a significant increase in the cumulative gas production. This phenomenon further shows that adsorbed gas constitutes an important source of shale gas that ensures a long period of stable production during the later stage of production.

Production modeling
Physical model. The wide range of pore sizes present in a shale reservoir leads to complex shale gas flow mechanisms that involve flow paths ranging from the molecular scale to the macroscale 38 . Shale gas production is a result of the combined action of multiple cross-scale migration mechanisms. This process can be divided into three main stages. During the first stage, the free gas in the fractures (including the secondary and artificial fractures) migrates toward the well shaft due to the decrease in the pressure in the nearby area. At the second stage, the free gas in the matrix flows toward the natural or artificial fractures due to the difference in the pressure between the matrix and fracture systems. During the third stage, the decrease in the pressure within the pores in a) Cumulative gas production versus production time (b) Cumulative gas production versus modified pressure Figure 10. Variation in the total, adsorbed, and free gas production. www.nature.com/scientificreports/ the matrix promotes the occurrence of microflow mechanisms (e.g., the diffusion and slip of the free gas and the diffusion of the adsorbed gas on the surface), resulting in an increased gas flowability in the matrix. When the pore pressure in the matrix decreases to the critical desorption pressure for the adsorbed gas, the adsorbed gas begins to contribute to gas production by desorbing as free gas into the pores in the matrix. Complex factors (e.g., natural fractures, brittle mineral content, and in situ stress) and various forms of hydraulic fractures induced by stimulation present challenges for accurately obtaining the characteristic parameters of a shale reservoir that describe its matrix and fracture heterogeneity. Therefore, it is assumed in this study that the forms of the main fractures in a reservoir can be characterized with a flat-plate double-wing model 39 . In addition, the fracture network and matrix region are treated as equivalent media. Some regions of the reservoir between and outside the main fractures are unfractured, resulting in a substantial interregional difference in the seepage behavior. Considering this factor, a quarter of a single main fracture is modeled. The model is divided into five regions: the main fracture in the internal region, the fracture network (region 1), the inter-main-fracture matrix (region 2), the unfractured matrix (region 3), and the matrix (region 4), as shown in Fig. 11.
The assumptions used to establish the model are detailed below: 1. Shale gas flows independently from the external regions (3 and 4) into the internal regions (1 and 2). The gas in region 2 converges toward the well shaft through region 1 and the main fracture. The migration of single-phase methane in different seepage fields is a 1D seepage process. Let y e denote the half-width of the shale reservoir and well length (marked as L) represents shale reservoir length in the model. The external boundaries are closed and homogeneous with a uniform thickness. 2. The adsorbed gas in the unfractured matrix (regions 2, 3, and 4) follows the supercritical Langmuir adsorption equation. The diffusion and transition flow of the free gas and the stress sensitivity of the matrix are considered. 3. The free gas in the equivalent fracture network migrates toward the main fracture by viscous flow, and stress sensitivity is nonnegligible. 4. The main fracture in the internal region is evenly distributed and vertically symmetrical with a uniform length and completely runs through the reservoir in the vertical direction. Let y F , w, L F , and d/2 denote the half-length and width of the main fracture, the spacing between main fracture clusters, and the half-width of a single hydraulic fracture, respectively. The gas flow behavior in the main fracture in the internal region follows Darcy's law. As the closure pressure increases, the proppant in the main fracture gradually undergoes changes (e.g., the proppant becomes elastically embedded into the wall of the main fracture and crushed), reducing the effective flow space in the main fracture. 5. The crossflow of gas from one seepage region to another is unsteady. Gas production is an isothermal seepage process. The effects of gravity and capillary forces are neglected. The bottomhole pressure in the gas well Figure 11. Simplified five-region physical model for a multicluster hydraulically fractured horizontal shale gas well. www.nature.com/scientificreports/ decreases from the original reservoir pressure to the predetermined output pressure during the initial stage of production, followed by constant-pressure production.

High-pressure physical property parameters
Because the high-pressure physical property parameters of the gas in the governing equation for seepage vary relatively considerably with temperature and pressure, the nonlinear effects cannot be neglected. A pseudopressure and a pseudotime are introduced in this study to linearize the governing equation to simplify the procedure used to solve it.
The pseudopressure is expressed as follows 38 : The pseudotime is expressed as follows 38 : where µ i is gas viscosity at initial formation pressure condition, Pa · s ; C ti is comprehensive compression factor at initial formation pressure condition, Pa −1 ; µ(p) is gas viscosity at formation pressure of p,Pa · s ; C ti is comprehensive compression factor at formation pressure of p ,, Pa −1 .

Gas density
The density of the free-phase gas in a seepage region can be derived from the state equation for a real gas, as shown below: where M is molar mass of methane (16 g/mol); R is universal gas constant (8.314 J/mol/K); T-reservoir temperature, K.

Supercritical adsorption-desorption model for gas
Analysis of the isothermal adsorption and development simulation test results presented earlier reveals the following phenomena that occur during the extraction of gas from a shale gas reservoir through depressurization. When the formation pressure is higher than the critical desorption pressure, adsorbed gas remains basically unextracted, while the gas well primarily produces free gas. When the formation pressure is lower than the critical desorption pressure, the gas originally in the adsorbed state desorbs and complements the free gas extracted from the reservoir, resulting in a decreased rate at which the pressure decreases and therefore ensuring a long period of stable production from the shale gas well during the later stage. There is a point on the high-pressure isothermal adsorption curve that corresponds to the maximum adsorption. When the pressure is higher than the critical pressure at this point, the adsorption is negatively correlated with the pressure 40,41 . The conventional Langmuir adsorption equation is unsuitable for describing high-pressure isothermal adsorption curves for shale reservoirs (Fig. 1). In this section, the following excess high-pressure isothermal adsorption model established based on the adsorbed-phase volume theory 40 is used: where q ad is supercritical excess amount of adsorbed shale gas per unit volume of the matrix, kg/m 3 ; Z sc is compressibility factor for an ideal gas, dimensionless.

Apparent permeability model
During the later stage of production of a shale gas reservoir, the flow regime in the micro-and nanopores in the matrix is not controlled by a single migration mechanism but instead is a result of the coupling of multiple flow mechanisms (e.g., viscous flow, diffusion, and transition flow). Establishing a multiscale, multiflow-regime apparent permeability model to reveal micro-and nanoscale gas flow patterns is a key scientific problem associated with production modeling. In this study, the stress sensitivity of the matrix and the desorption of gas in the matrix are coupled. Then, the apparent permeability model established by Wu et al. 42 for micro-and nanopores in shale matrices based on molecular dynamics theories is used to superpose the viscous flow and Knudsen diffusion of real gas molecules with weighting coefficients: www.nature.com/scientificreports/ The effective Knudsen number is given below: The effective hydraulic flow radius that accounts for stress sensitivity and desorption is expressed as follows: The coverage of gas adsorbed onto the wall surface of the pores in the matrix is calculated below: where K ma is effective apparent matrix permeability, m 2 ; K ne is effective Knudsen number, dimensionless;K m0 is matrix intrinsic permeability, m 2 ; C g is gas compressibility coefficient, Pa −1 ; D is gas diffusion coefficient, m 2 /s; C d is compressibility factor for desorbed gas in the matrix, Pa -1 ; C m is matrix compressibility coefficient, Pa −1 ; D is gas diffusion coefficient, m 2 /s; C d is compressibility factor for desorbed gas in the matrix, Pa -1 ; is length of the free path of methane molecules at any arbitrary temperature and pressure, m; r e is effective hydraulic flow radius, m; r 0 is initial hydraulic flow radius, m; p m is shale matrix pressure, Pa; γ m is shale matrix stress sensitivity, Pa −1 ; p is average pressure in the reservoir, Pa; K B is Boltzmann constant (1.38065 × 10 -23 J/K); d is collision diameter of methane molecules, m; d CH 4 is diameter of methane molecules, m.

Pressure sensitivity effect
When gas is extracted from a fracture network, the non-proppant-supported secondary fractures may be closed due to stress sensitivity, resulting in a significantly lower permeability of the secondary fracture network. Therefore, an empirical stress sensitivity model in an exponential form 41 , given below, is used in this study to characterize the effects of the stress sensitivity of secondary fractures on gas flowability.
where K f is permeability that accounts for pressure sensitivity, m 2 ; K fi is permeability at the initial time, m 2 ; γpseudo permeability modulus, Pa 2 /(Pa 2 ·s).
Based on the dimensionless parameters defined in Table 2, the system of governing equations for seepage is converted to dimensionless seepage equations to solve the dimensionless productivity model.
Where ψ i is initial formation pseudo-pressure, Pa 2 /(Pa · s); ψ is the formation pseudo-pressure, Pa 2 /(Pa · s); K F is hydraulic fracture permeability, m 2 ; ϕ 2 is porosity of the matrix zone 2, dimensionless; ϕ 1 is porosity of fracture network, dimensionless; ϕ 3 is porosity of matrix zone 3. m 2 ; C t3 is comprehensive compressibility of matrix zone 3, Pa −1 ; C t2 is comprehensive compressibility of matrix zone 2, Pa −1 ; C t1 is comprehensive compressibility of fracture network, Pa −1 ; t a is pseudo-time, s; L f is hydraulic fracture spacing, m; A cw is wellbore crossflow area, m 2 ; ϕ F is porosity of hydraulic fracture, m; C tF is comprehensive compressibility of hydraulic fracture, Pa −1 ; K 3a is apparent permeability of matrix zone 3, m 2 ; K 1 is fracture network permeability, m 2 ; K 2 is matrix zone 2 permeability, m 2 ;d is the width of stimulated reservoir zone, m.
Mathematical model. 1. Governing equation for seepage in the matrix in region 4 www.nature.com/scientificreports/ The gas in the matrix in region 4 converges along the y-direction into the matrix in region 2 by unsteady crossflow. Considering the supercritical desorption of the adsorbed gas and the Knudsen diffusion and viscous flow of the gas in the nanopores in the shale matrix, a closed external boundary and a continuous pressure at the internal boundary are used 43 .
Based on the law of conservation of mass, a governing equation for seepage in the matrix is obtained, as shown below: The seepage velocity of the gas in the matrix can be derived, as shown below: Adding the excess gas adsorption in the form of the compressibility factor for gas in the desorbed state into the compressibility factor for gas in the free state yields a comprehensive compressibility factor for gas in the matrix in the external region, C t4 : Combining Eqs. (19), (20), and (21) and substituting the result into Eq. (18) transforms it to Substitution of Eqs. (17) and (22) into Eq. (16) gives where ϕ 4 is porosity of matrix zone 4, dimensionless; C t4 is the matrix comprehensive compressibility coefficient, Pa −1 ; P 4 is pore pressure in matrix zone 4, Pa; t is production duration, s; K 4a is apparent permeability of matrix zone 4, m 2 ; µ 4 is gas viscosity, Pa · s ; C g4 is gas compressibility, Pa −1 ; C d4 is modified supercritical desorption gas compression coefficient of matrix zone 4, Pa −1 ; C f4 is formation compressibility, Pa −1 ; ϕ 4 is porosity of the matrix zone 4, dimensionless; K 4a is apparent permeability of matrix zone 4, m 2 ; K 4i is intrinsic permeability of matrix zone 4, m 2 .
For simplicity, Eq. (23) is transformed into a pseudo-function with dimensionless variables, which is subsequently simplified to yield Eq. (24): The pressure in the seepage field in the external region at the initial time is equivalent to the initial formation pressure, as shown below: Equations (25)- (27) are transformed into the following forms of dimensionless boundary conditions: Similar to the gas seepage pattern in the matrix in region 4, gas in the matrix in region 3 converges along the y-direction toward the matrix in region 1 by Knudsen diffusion, viscous flow, and unsteady crossflow. A dimensionless governing equation for seepage in the matrix in region 3 can be derived from the mass conservation, motion, and state equations, as shown below: Because the external boundary is closed and a continuous pressure is present at the internal boundary, the dimensionless initial and boundary conditions expressed in Eqs. (40)-(42) hold: www.nature.com/scientificreports/ A dimensionless pseudopressure solution in the Laplace form for the matrix in region 3, shown in Eq. (43), is obtained through a Laplace transform using a procedure similar to that used to solve the dimensionless governing equation for flow in the matrix in region 4.

Governing equation for seepage in the matrix in region 2
Gas in the matrix in region 2 flows into the matrix in region 1 along the x-direction by diffusion and viscous flow. In addition, the external boundary is closed, while a continuous pressure is present at the internal boundary in regions 2 and 1. Considering these factors, a dimensionless governing equation for seepage in the matrix in region 2 is established, as shown below: The external boundary allows no flow. The pseudopressure in the fracture network at the internal boundary is identical to that in the main fracture. The pseudopressure in the seepage field in the fracture network at the initial time is equivalent to the initial formation pseudopressure. Thus, the initial boundary conditions can be expressed as follows: A Laplace transform is performed to derive a dimensionless pseudopressure solution for the governing equation for flow in the matrix in region 2:  Considering that the gas well produces gas at a constant bottomhole pressure and that the external boundary is closed, the following dimensionless initial boundary conditions hold: A pseudopressure solution in the Laplace space for the derived dimensionless governing equation for seepage in the main fracture is given below: where Let N denote the number of hydraulic fractures in the shale gas well. Then, the dimensionless gas production rate at the bottom of the shale gas well can be expressed in the Laplace space as follows: Subsequently, a real-space semianalytical dimensionless gas production rate is obtained through a Stehfest numerical inversion 40 . Finally, Newton's iteration method can be used to determine the real-space gas production from the gas well at a constant pressure.

Model validation and prediction
Fitting physical simulation test results. To examine the effectiveness of the model, the matrix part of the model is validated based on the physical simulation test results. The developed mathematical model is first simplified by assuming that the hydraulic fracture and fracture network regions have infinite flow conductivity and by neglecting the gas flow in the reservoir in the y-direction in the model. Then, a 1D gas production pattern is simulated for the matrix at the core scale. Subsequently, the gas production test data for the matrix core are fitted. Finally, the model-fitted parameter (Table 3) is inverted to validate the model. (* indicates a model-fitted parameter). www.nature.com/scientificreports/ Figures 12 and 13 compare the measured and model-fitted gas production curves. Overall, a relatively high goodness-of-fit can be observed between the gas production curves yielded by the developed five-region productivity model and the test. The model yields a cumulative gas production of 4756 mL, which differs by less than 5% from the actual value for the gas well after 3000 days of production. This finding validates the soundness and feasibility of the coupling of the high-pressure adsorption, apparent permeability, and stress sensitivity models in the developed model, while showing that the model can be used to theoretically analyze the productivity of shale gas wells and predict their production patterns. Note that there is a marked difference between the actual and model-fitted production curves for a period of time during the initial stage, which can be primarily ascribed to the overly fast decrease in the outlet pressure during production. A relatively high goodness-of-fit can be observed between the measured and model-fitted curves for the stable period of production during the later stage.
Single-well simulation and production prediction. The productivity model developed and validated based on the physical simulation test results earlier was used to fit the dynamic production data for a multistagefractured horizontal well (NH3-5-1) in the Changning Shale Gas Demonstration Area in the Sichuan Basin and to analyze its full-life-cycle dynamic production. Table 4 summarizes the relevant geological parameters and well completion parameters. (* indicates a model-fitted parameter).
Observation of the model-fitted gas production curves in Figs. 14 and 15 reveals the following. The production of the gas well lasted for 1500 days. The daily production fluctuated relatively appreciably but, overall, decreased   www.nature.com/scientificreports/ in a relatively markedly monotonic pattern. Over the first 900 days, the daily gas production rate decreased from 300,000 to approximately 50,000 m 3 , while the cumulative gas production reached 0.98 × 10 8 m 3 . The gas well had basically entered the stable production stage during this time period. Generally similar trends can be observed between the gas production curves yielded by the developed productivity model and the actual dynamic production curves. At first months of production, production data is lower than the model result, which can be possibly attributable to the fracturing fluid back flow, excessive choking and/or lack of takeaway capacity [44] . The model yields a cumulative gas production of 0.98 × 10 8 m 3 , which differs by approximately 1.27% from the actual value after 1500 days of production. These findings show that the developed productivity model is reliable and can be used to accurately generate productivity and dynamic production pattern estimates for gas wells.  Figure 14. Fit for the daily gas production of well NH3-5-1. www.nature.com/scientificreports/ Figures 16 and 17 show that when the non-linear gas flow mechanism is not considered in the production model, the predicted 20-year-production predicted by the model can be underestimated by 11.3%, indicating diffusion and desorption is significant to the later-stage-productivity; the calculated EUR of considering stress sensitivity is 7.5% lower than that model predictions without involving stress sensitivity. Likewise, the model considering the high-pressure adsorption predicts that EUR is 1.87% higher than that of only considering Langmuir adsorption. If the above comprehensive factors are considered in the model, the EUR is about 2.5% higher than that of the linear production model. To sum up, the comprehensive factors considered in the production model is conducive to enhance the accuracy of predicting the medium and long-term productivity. It aims to provide reliable reference to accurately estimate shale gas well production at different production durations. The developed productivity model was used to investigate well NH3-5-1 in terms of development dynamics and the variation in adsorbed and free gas production over a 20-year-period. Figure 18 shows the results. For   www.nature.com/scientificreports/ the initial stage of production (i.e., approximately the first two years of production), the cumulative gas production curve basically overlaps the free gas production curve, while the adsorbed gas production is relatively low. Therefore, free gas accounts for the majority of the gas produced during the initial stage. For the middle and later stages of production, the cumulative gas production curve gradually deviates from the free gas production curve, while the adsorbed gas production increases appreciably. The EUR of the gas well is 156 million m 3 , of which adsorbed and free gas account for 34 and 122 million m 3 , respectively. Our observation of Fig. 19 shows that the proportion of adsorbed gas increases gradually as production proceeds, which differs completely from the pattern of variation in the proportion of free gas. Ultimately, adsorbed gas accounts for as much as 20% of the total gas produced from the well. Analysis of the variation in the gas production yielded by the model with the formation pressure reveals the following. At an average formation pressure greater than 35 MPa, the cumulative gas production is basically linearly related with the apparent formation pressure. The slope of the corresponding straight line is 1.81 million m 3 /MPa. Under this condition, the cumulative gas production curve also basically coincides with the free gas production curve, suggesting that the gas produced from the well is basically in the free state. As the formation pressure decreases below 35 MPa, the adsorbed gas begins to be extracted and contributes more to gas production. Consequently, the cumulative gas production curve begins to deviate from the straight line and gradually bends upward. As the apparent formation pressure decreases, the extent to which each of the adsorbed and cumulative gas production curves bends upwards increases. The average slopes of these two curves are 3.6 and 4.02 million m 3 /MPa. The adsorbed gas content of the gas produced from the well plays a relatively crucial role. At the final stage of production (i.e., at a reservoir pressure of approximately 10 MPa), the cumulative adsorbed gas production reaches 340 million m 3 , accounting for approximately 21.79% of the EUR. Therefore, the role of adsorbed gas in the dynamic production process becomes increasingly prominent during the later stage when the reservoir pressure is relatively low. The productivity of the reservoir relies predominantly on the desorption of adsorbed gas during the later stage. Observation of Figs. 20 and 21 shows the following. During the first year of production of the gas well, the pressure in the reservoir averages approximately 40 MPa. Free gas accounts for the majority of the produced gas, while adsorbed gas remains basically unextracted and contributes only 4.1% of the total annual gas production. In addition, 40.8% of the EUR is recovered during the first year. Over time, the pressure in the reservoir decreases, and the annual gas production decreases in a notably L-shaped pattern. Free gas is extracted in large quantities, while the adsorbed gas production increases gradually. In the seventh year, adsorbed gas accounts for as much as 50% of the produced gas, while the pressure in the reservoir decreases below 20 MPa, with more than 80% of the EUR having been recovered. During the middle and later stages, annual gas production remains at a constant low rate of 100 million m 3 , the majority of which is supplied by the extraction of adsorbed gas. In Figure 19. Variation in the proportions of adsorbed and free gas with production time. Figure 20. Variation in the annual gas production calculated for the well with production time. www.nature.com/scientificreports/ the 20th year, the remaining formation pressure averages 10 MPa, an approximately 85% decrease compared with its initial level. The gas production is calculated to be 155.6 million m 3 . The cumulative amount of adsorbed gas produced from the well makes up 21% of the EUR. Therefore, the potential of adsorbed gas should be fully exploited during actual production to improve daily and cumulative gas production and ensure a long period of stable production from the gas well.

Discussions
This paper designed a full-diameter shale core development physical simulation under high pressure conditions for over 1000 days, taking place of the small shale core short term extraction experiments in previous studies, which makes it reliable to simulate field-scale shale gas production. However, high temperature condition was not included in the proposed experiment in the paper, which is a future research direction in physical simulation of high pressure-high temperature gas extraction for shale gas reservoirs.
In addition, the proposed model is incorporated in comprehensive gas nonlinear effects, including its physical properties at high pressures, supercritical desorption characteristics, and multiple flow mechanisms. This model is more convincing in predicting and analyzing shale gas production decline trend in different production duration and produced gas compositions. For instance, when the non-linear gas flow mechanism is not considered in the production model, the predicted 20-year-production can be underestimated by 11.3%; the calculated EUR of considering stress sensitivity is 7.5% lower than that model predictions without involving stress sensitivity. Likewise, the model considering the high-pressure adsorption predicts that EUR is 1.87% higher than that of only considering Langmuir adsorption. If the above comprehensive factors are considered in the model, the EUR is about 2.5% higher than that of the linear production model. For that reason, the more comprehensive factors considered in the production model is conducive to enhance the accuracy of predicting the medium and long-term productivity. It aims to provide reliable reference to accurately estimate shale gas well production at different production durations.
However, the proposed production model is related to only describe single gas flow process, without considering the gas-water two phases transport in the backflow period. The gas-water flow characterization is exactly to be involved in the production model in the future works.
To sum up, the findings from study can help us explicitly understand and predict shale gas process out of subsurface formations under high geological pressure conditions, dominated by multiple flow mechanisms.

Conclusions
In this study, a 3343-day ultralong-production-cycle development simulation was conducted with methane as the test gas. As indicated by the pressure and gas production rate curve, the obtained gas accumulation and production characteristics were highly consistent with those observations in the reservoir. During gas-well production, the low gas flowability in the shale matrix led to relatively slow pressure propagation. As a result, only 71% of the EUR was recovered after 3343 days. The production decreased rapidly during the initial stage, followed by a long period of low and stable production during the later stage. This phenomenon was basically consistent with that inreal gas well. A critical desorption pressure (approximately 12 MPa) was found for the extraction of adsorbed gas. Adsorbed gas did not contribute to gas production until the pressure decreased below the critical desorption pressure.
A five-region seepage model that accounts for high-pressure supercritical adsorption, apparent permeability, and stress sensitivity was developed and subsequently validated based on the physical simulation test results. Then it was proved that more gas flow nonlinear effects considered in the proposed model can effectively promote higher prediction accuracy in the medium and long-term productivity of the gas well.
In addition, the production decline analysis were conducted on the specifically single well production data with the production model, which illustrated that free gas accounts for the majority of the produced gas during the initial extraction stage. In the first year, 40.8% of the EUR is recovered, with free gas accounting for more than 90% of the produced gas. During the later production stage, the adsorbed gas is a primary source of gas. www.nature.com/scientificreports/ Adsorbed gas contributes as much as 50% of the total gas produced in the seventh year. The 20-year-cumulative adsorbed gas production makes up 21% of the EUR. Adsorbed gas is a major factor that can ensure a long period of stable production during the later stage. During the later stage of production, adsorbed gas begins to desorb, resulting in expanded pore flow channels and an increase in gas flowability. The decreased pressure also enhances gas diffusion and thereby correspondingly increases gas flowability. Therefore, shale gas should be extracted from a reservoir by means of depressurization. The bottomhole pressure should be decreased below the critical desorption pressure as fast as possible to enable effective extraction of the adsorbed gas. However, in terms of the gas-well development, gas should be produced at a controlled pressure for a certain period of time during the initial stage to reduce the negative impact of fracturing fluid flowback on the production, fracture and fracture network closure, and sand production. Subsequently, the pressure should be reduced as fast as possible to exploit the adsorbed gas, aimed at improving daily gas production and ultimate recovery.